Haptic aesthetics and bodily properties of Ori Gersht’s digital art a behavioral and eye tracking study

1 Latency of First Fixation ———————–

Normalized Latency of First Fixtion Database Import

FixLatency <- readr::read_delim("FixLatency.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)

Parsed with column specification:
cols(
  .default = col_double(),
  RecMedia = col_character(),
  Cross = col_character(),
  Question = col_character(),
  Resolution = col_character(),
  Magnitude = col_character(),
  RecordingNameII = col_character(),
  MediaName = col_character(),
  LeftHit = col_character(),
  RightHit = col_character(),
  AOI = col_character()
)
See spec(...) for full column specifications.

FixLatency <- subset(FixLatency, !(FixLatency$Question == "Baseline"))
FixLatency$AOI <- as.factor(FixLatency$AOI)
FixLatency$Magnitude <- as.factor(FixLatency$Magnitude)
levels(FixLatency$Magnitude) <- c("Detail", "Entire")
print(FixLatency)

Latency Definitive Model

Lat <- lme4::lmer(LatencyNorm ~  AOI + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixLatency)

Latency Overall Statistics

AIC(Lat)

[1] 28819.81

MuMIn::r.squaredGLMM(Lat)

'r.squaredGLMM' now calculates a revised statistic. See the help page.

            R2m        R2c
[1,] 0.02223695 0.09513032

Latency Summary & ANOVA

summary(Lat)

Linear mixed model fit by REML ['lmerMod']
Formula: LatencyNorm ~ AOI + AOI:Magnitude + (1 | AOI:RecordingNameII)
   Data: FixLatency

REML criterion at convergence: 28807.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.3826 -0.4953 -0.2142  0.1213 10.4241 

Random effects:
 Groups              Name        Variance Std.Dev.
 AOI:RecordingNameII (Intercept) 0.106    0.3256  
 Residual                        1.316    1.1471  
Number of obs: 9203, groups:  AOI:RecordingNameII, 58

Fixed effects:
                         Estimate Std. Error t value
(Intercept)               0.83815    0.06500  12.894
AOIRight                  0.21382    0.09193   2.326
AOILeft:MagnitudeEntire  -0.08450    0.03376  -2.503
AOIRight:MagnitudeEntire  0.16059    0.03388   4.740

Correlation of Fixed Effects:
            (Intr) AOIRgh AOIL:M
AOIRight    -0.707              
AOILft:MgnE -0.260  0.184       
AOIRght:MgE  0.000 -0.183  0.000

aovLat <- aov(LatencyNorm ~ AOI + AOI:Magnitude + Error(RecordingNameII/(Resolution + Magnitude)), data = FixLatency)
summary(aovLat)

Error: RecordingNameII
              Df Sum Sq Mean Sq F value Pr(>F)  
AOI            1   5.94   5.941   1.414 0.2455  
AOI:Magnitude  2  30.25  15.123   3.600 0.0423 *
Residuals     25 105.03   4.201                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: RecordingNameII:Resolution
              Df Sum Sq Mean Sq F value Pr(>F)
AOI            1   0.62  0.6222   0.501  0.485
AOI:Magnitude  2   0.34  0.1685   0.136  0.874
Residuals     26  32.27  1.2411               

Error: RecordingNameII:Magnitude
              Df Sum Sq Mean Sq F value Pr(>F)
AOI            1   2.39   2.392   1.717  0.202
AOI:Magnitude  2   5.01   2.507   1.799  0.185
Residuals     26  36.23   1.393               

Error: Within
                Df Sum Sq Mean Sq F value   Pr(>F)    
AOI              1    257  257.29  182.81  < 2e-16 ***
AOI:Magnitude    1     34   33.81   24.02 9.69e-07 ***
Residuals     9114  12828    1.41                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Latency Robustlmm

require(lme4)
require(robustlmm)
robustLat <- robustlmm::rlmer(LatencyNorm ~  AOI + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixLatency, rho.sigma.e = psi2propII(smoothPsi, k = 2.28), rho.sigma.b = psi2propII(smoothPsi,k = 2.28))

summary(robustLat)

Robust linear mixed model fit by DAStau 
Formula: LatencyNorm ~ AOI + AOI:Magnitude + (1 | AOI:RecordingNameII) 
   Data: FixLatency 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.7333 -0.5259 -0.1955  0.4628 18.2600 

Random effects:
 Groups              Name        Variance Std.Dev.
 AOI:RecordingNameII (Intercept) 0.04645  0.2155  
 Residual                        0.45694  0.6760  
Number of obs: 9203, groups: AOI:RecordingNameII, 58

Fixed effects:
                         Estimate Std. Error t value
(Intercept)               0.63699    0.04351  14.640
AOIRight                  0.20903    0.06154   3.397
AOILeft:MagnitudeEntire  -0.03172    0.02041  -1.554
AOIRight:MagnitudeEntire  0.07446    0.02048   3.636

Correlation of Fixed Effects:
            (Intr) AOIRgh AOIL:M
AOIRight    -0.707              
AOILft:MgnE -0.235  0.166       
AOIRght:MgE  0.000 -0.166  0.000

Robustness weights for the residuals: 
 7695 weights are ~= 1. The remaining 1508 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0737  0.3440  0.5920  0.5920  0.8520  0.9990 

Robustness weights for the random effects: 
 46 weights are ~= 1. The remaining 12 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.582   0.726   0.836   0.825   0.961   0.996 

Rho functions used for fitting:
  Residuals:
    eff: smoothed Huber (k = 1.345, s = 10) 
    sig: smoothed Huber, Proposal II (k = 2.28, s = 10) 
  Random Effects, variance component 1 (AOI:RecordingNameII):
    eff: smoothed Huber (k = 1.345, s = 10) 
    vcp: smoothed Huber, Proposal II (k = 2.28, s = 10) 

Comparing lmer-rlmer models

campareLat <- robustlmm::compare(Lat, robustLat, show.rho.functions = FALSE)

Get the bootstrapped confidence intervals for the model parameters

ciLat <- lme4::confint.merMod(Lat, method = "boot", boot.type = "perc", nsim = 10000, oldNames = F)

ciLat

                                         2.5 %      97.5 %
sd_(Intercept)|AOI:RecordingNameII  0.26053283  0.38945797
sigma                               1.13069224  1.16372164
(Intercept)                         0.70949223  0.96761043
AOIRight                            0.03295627  0.39482166
AOILeft:MagnitudeEntire            -0.15037329 -0.01928073
AOIRight:MagnitudeEntire            0.09381820  0.22762033

Latency Post Hoc

emmeans::emm_options(pbkrtest.limit = 9203)
emmLat <- emmeans::emmeans(Lat, pairwise ~ AOI:Magnitude, adjust="tukey")

emmLat

$emmeans
 Magnitude AOI   emmean     SE   df lower.CL upper.CL
 Detail    Left   0.838 0.0650 64.4    0.708    0.968
 Entire    Left   0.754 0.0650 64.4    0.624    0.883
 Detail    Right  1.052 0.0650 64.4    0.922    1.182
 Entire    Right  1.213 0.0651 64.6    1.083    1.342

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast                    estimate     SE     df t.ratio p.value
 Detail,Left - Entire,Left     0.0845 0.0338 9143.1  2.503  0.0595 
 Detail,Left - Detail,Right   -0.2138 0.0919   64.4 -2.326  0.1028 
 Detail,Left - Entire,Right   -0.3744 0.0920   64.5 -4.071  0.0007 
 Entire,Left - Detail,Right   -0.2983 0.0919   64.4 -3.245  0.0098 
 Entire,Left - Entire,Right   -0.4589 0.0920   64.5 -4.990  <.0001 
 Detail,Right - Entire,Right  -0.1606 0.0339 9143.1 -4.740  <.0001 

P value adjustment: tukey method for comparing a family of 4 estimates 

1.1 Latency of First Fixation Hierarchical Steps —-

Latency Nested Models

null <- nlme::gls(LatencyNorm ~ 1, data = FixLatency)
AOI <- lm(LatencyNorm ~ AOI, data = FixLatency)
AOI.intAOIM <- lm(LatencyNorm ~ AOI + AOI:Magnitude, data = FixLatency)
AOI.intAOIM.rand.intAOIM <- lme4::lmer(LatencyNorm ~ AOI + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixLatency)

Latency Nested Models AIC

MuMIn::AICc(null)

[1] 29542.21

MuMIn::AICc(AOI)

[1] 29358.97

MuMIn::AICc(AOI.intAOIM)

[1] 29336.98

MuMIn::AICc(AOI.intAOIM.rand.intAOIM)

[1] 28819.82

Latency Nested Models BIC

BIC(null)

[1] 29556.46

BIC(AOI)

[1] 29380.35

BIC(AOI.intAOIM)

[1] 29372.61

BIC(AOI.intAOIM.rand.intAOIM)

[1] 28862.58

Latency Nested Models R squared

MuMIn::r.squaredLR(null)

[1] 2.775558e-15
attr(,"adj.r.squared")
[1] 2.892417e-15

MuMIn::r.squaredGLMM(AOI)

            R2m        R2c
[1,] 0.01918762 0.01918762

MuMIn::r.squaredGLMM(AOI.intAOIM)

           R2m       R2c
[1,] 0.0219494 0.0219494

MuMIn::r.squaredGLMM(AOI.intAOIM.rand.intAOIM)

            R2m        R2c
[1,] 0.02223695 0.09513032

Latency Nested Models Likelihood ratio test

lmtest::lrtest(null, AOI, AOI.intAOIM, AOI.intAOIM.rand.intAOIM)

Likelihood ratio test

Model 1: LatencyNorm ~ 1
Model 2: LatencyNorm ~ AOI
Model 3: LatencyNorm ~ AOI + AOI:Magnitude
Model 4: LatencyNorm ~ AOI + AOI:Magnitude + (1 | AOI:RecordingNameII)
  #Df LogLik Df   Chisq Pr(>Chisq)    
1   2 -14769                          
2   3 -14676  1 185.237  < 2.2e-16 ***
3   5 -14664  2  25.997  2.264e-06 ***
4   6 -14404  1 519.159  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

1.2 Latency of First Fixation Plots —————–

Interaction AOI * MAgnitude

Create Database with Latency fitted values (representing the effect on the scale of the linear predictor)

efLat <- effects::effect("AOI:Magnitude", Lat)
summary(efLat)

 AOI*Magnitude effect
       Magnitude
AOI        Detail    Entire
  Left  0.8381502 0.7536477
  Right 1.0519745 1.2125632

 Lower 95 Percent Confidence Limits
       Magnitude
AOI        Detail    Entire
  Left  0.7107345 0.6262365
  Right 0.9245440 1.0850473

 Upper 95 Percent Confidence Limits
       Magnitude
AOI        Detail    Entire
  Left  0.9655658 0.8810589
  Right 1.1794050 1.3400791

efLat <- as.data.frame(efLat)

Create Database with Latency means by Subject

require(dplyr)

Carico il pacchetto richiesto: dplyr

Attaching package: <U+393C><U+3E31>dplyr<U+393C><U+3E32>

The following objects are masked from <U+393C><U+3E31>package:stats<U+393C><U+3E32>:

    filter, lag

The following objects are masked from <U+393C><U+3E31>package:base<U+393C><U+3E32>:

    intersect, setdiff, setequal, union

LatmeansSubj <- FixLatency %>% 
  group_by(AOI, Magnitude, RecordingNameII) %>% 
  summarise(RespM = mean(LatencyNorm), Respse =  sjstats::se(LatencyNorm))
print(LatmeansSubj)

1.2.0.1 Plot Latency fitted values

require(ggplot2)

Carico il pacchetto richiesto: ggplot2

Latbar <- ggplot(efLat, aes(AOI, fit, fill= Magnitude)) + 
  geom_point() + 
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4) + theme_bw(base_size=12)
Latbar <- ggplot(efLat, aes(AOI, fit, fill=Magnitude)) + geom_bar(stat="identity", position="dodge") +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4, position=position_dodge(width=0.9)) + 
  ylim(0,1.3) +
  labs(x= "AOI", y= "Normalized First Fixation Latency") +
  scale_fill_grey(name = "Magnitude") +
  ggtitle("AOI*Magnitude") +
  theme_classic() +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(Latbar)

1.2.0.2 Plot Latency Means by Subject

require(ggplot2)
LatsplitSUBJ <- ggplot(LatmeansSubj, aes(AOI, RespM, fill= Magnitude)) + 
    geom_bar(stat="identity", position="dodge") +
    facet_wrap(~ RecordingNameII) + 
    xlab("AOI") + 
    ylab("Normalized First Fixation Latency") +
    scale_fill_grey(name = "Magnitude")
print(LatsplitSUBJ)

1.2.0.3 Boxplot Latency

Latbox <- ggplot(data = FixLatency, aes(x = Magnitude, y = LatencyNorm)) +
  facet_wrap(. ~ AOI) + 
  geom_jitter(alpha = 0.1, color = "black") +
  geom_boxplot(alpha = 0) +
  theme_light() +
  labs(tag = "A", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized First Fixation Latency") +
  theme(plot.title = element_text(hjust = 0.5))
print(Latbox)

1.2.0.4 Violinplot Latency

Function: plot split violins

GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
                             
                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })
geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}

Violinplot Latency

Latsplitviolin <- ggplot(FixLatency, aes(x =Magnitude, y=LatencyNorm, fill=AOI)) + 
  geom_split_violin() +
  scale_fill_grey() +
  theme_light() +
  labs(tag = "A", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized First Fixation Latency") +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(Latsplitviolin)

2 Total Number of Fixations ———————-

Normalized Total Number of Fixations Database Import

FixNumber <- readr::read_delim("FixNumber.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)

Parsed with column specification:
cols(
  RecMedia = col_character(),
  RecordingNameII = col_character(),
  MediaName = col_character(),
  ImageStartTimestamp = col_double(),
  Cross = col_character(),
  Question = col_character(),
  Resolution = col_character(),
  Magnitude = col_character(),
  AOI = col_character(),
  NFixAOIQRM = col_double(),
  NFixAOIQRMNorm = col_double()
)

FixNumber <- subset(FixNumber, !(FixNumber$Question == "Baseline"))
FixNumber$AOI <- as.factor(FixNumber$AOI)
FixNumber$Magnitude <- as.factor(FixNumber$Magnitude)
levels(FixNumber$Magnitude) <- c("Detail", "Entire")
FixNumber$Resolution <- as.factor(FixNumber$Resolution)
print(FixNumber)

NFix Definitive Model

NFix <- lme4::lmer(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixNumber)

NFix Overall Statistics

AIC(NFix)

[1] 9274.25

MuMIn::r.squaredGLMM(NFix)

            R2m       R2c
[1,] 0.03773117 0.2759189

NFix Summary & ANOVA

summary(NFix)

Linear mixed model fit by REML ['lmerMod']
Formula: NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude +      (1 | AOI:RecordingNameII)
   Data: FixNumber

REML criterion at convergence: 9260.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4775 -0.6892 -0.0227  0.6324  4.4338 

Random effects:
 Groups              Name        Variance Std.Dev.
 AOI:RecordingNameII (Intercept) 0.05083  0.2254  
 Residual                        0.15451  0.3931  
Number of obs: 9278, groups:  AOI:RecordingNameII, 58

Fixed effects:
                          Estimate Std. Error t value
(Intercept)               1.049192   0.042848  24.487
AOIRight                 -0.038203   0.060320  -0.633
ResolutionLow            -0.031963   0.008162  -3.916
MagnitudeEntire           0.078235   0.011543   6.778
AOIRight:MagnitudeEntire -0.206100   0.016324 -12.626

Correlation of Fixed Effects:
            (Intr) AOIRgh RsltnL MgntdE
AOIRight    -0.704                     
ResolutinLw -0.095  0.000              
MagnitdEntr -0.135  0.096  0.000       
AOIRght:MgE  0.095 -0.135  0.000 -0.707

aovNFix <- aov(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + Error(RecordingNameII/(Resolution + Magnitude)), data = FixNumber)
summary(aovNFix)

Error: RecordingNameII
           Df Sum Sq Mean Sq F value Pr(>F)
Resolution  1  0.109  0.1091   0.178  0.676
Residuals  27 16.518  0.6118               

Error: RecordingNameII:Resolution
           Df Sum Sq Mean Sq F value   Pr(>F)    
Resolution  1  2.369  2.3690  15.802 0.000473 ***
Magnitude   1  0.004  0.0037   0.025 0.875849    
Residuals  27  4.048  0.1499                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: RecordingNameII:Magnitude
          Df Sum Sq Mean Sq F value Pr(>F)   
Magnitude  1  1.428  1.4282   12.76 0.0013 **
Residuals 28  3.133  0.1119                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
                Df Sum Sq Mean Sq F value Pr(>F)    
AOI              1   46.3   46.29   228.2 <2e-16 ***
AOI:Magnitude    1   24.7   24.65   121.5 <2e-16 ***
Residuals     9189 1864.4    0.20                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

NFix Robustlmm

require(lme4)
require(robustlmm)
robustNFix <- robustlmm::rlmer(NFixAOIQRMNorm ~  AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixNumber, rho.sigma.e = psi2propII(smoothPsi, k = 2.28), rho.sigma.b = psi2propII(smoothPsi, k = 2.28))

summary(robustNFix)

Robust linear mixed model fit by DAStau 
Formula: NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude +      (1 | AOI:RecordingNameII) 
   Data: FixNumber 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5441 -0.6613 -0.0117  0.6541  4.5708 

Random effects:
 Groups              Name        Variance Std.Dev.
 AOI:RecordingNameII (Intercept) 0.05763  0.2401  
 Residual                        0.15106  0.3887  
Number of obs: 9278, groups: AOI:RecordingNameII, 58

Fixed effects:
                          Estimate Std. Error t value
(Intercept)               1.038044   0.046648  22.253
AOIRight                 -0.034469   0.065710  -0.525
ResolutionLow            -0.029371   0.008277  -3.549
MagnitudeEntire           0.082481   0.011705   7.046
AOIRight:MagnitudeEntire -0.222657   0.016554 -13.451

Correlation of Fixed Effects:
            (Intr) AOIRgh RsltnL MgntdE
AOIRight    -0.704                     
ResolutinLw -0.089  0.000              
MagnitdEntr -0.125  0.089  0.000       
AOIRght:MgE  0.089 -0.126  0.000 -0.707

Robustness weights for the residuals: 
 7385 weights are ~= 1. The remaining 1893 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.294   0.672   0.823   0.792   0.939   0.999 

Robustness weights for the random effects: 
 47 weights are ~= 1. The remaining 11 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.590   0.692   0.809   0.807   0.960   0.974 

Rho functions used for fitting:
  Residuals:
    eff: smoothed Huber (k = 1.345, s = 10) 
    sig: smoothed Huber, Proposal II (k = 2.28, s = 10) 
  Random Effects, variance component 1 (AOI:RecordingNameII):
    eff: smoothed Huber (k = 1.345, s = 10) 
    vcp: smoothed Huber, Proposal II (k = 2.28, s = 10) 

Comparing lmer-rlmer models

campareNFix <- robustlmm::compare(NFix, robustNFix, show.rho.functions = FALSE)

Get the bootstrapped confidence intervals for the model parameters

ciNFix <- lme4::confint.merMod(NFix, method = "boot", boot.type = "perc", nsim = 10000, oldNames = F)

ciNFix

                                         2.5 %      97.5 %
sd_(Intercept)|AOI:RecordingNameII  0.26053283  0.38945797
sigma                               1.13069224  1.16372164
(Intercept)                         0.70949223  0.96761043
AOIRight                            0.03295627  0.39482166
AOILeft:MagnitudeEntire            -0.15037329 -0.01928073
AOIRight:MagnitudeEntire            0.09381820  0.22762033

NFix Post Hoc

emmeans::emm_options(pbkrtest.limit = 9278)
emmNFix <- emmeans::emmeans(NFix, pairwise ~ AOI:Magnitude, adjust="tukey") 

emmNFix

$emmeans
 AOI   Magnitude emmean     SE   df lower.CL upper.CL
 Left  Detail     1.033 0.0427 58.1    0.948    1.119
 Right Detail     0.995 0.0427 58.1    0.910    1.080
 Left  Entire     1.111 0.0427 58.1    1.026    1.197
 Right Entire     0.867 0.0427 58.1    0.782    0.953

Results are averaged over the levels of: Resolution 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast                    estimate     SE     df t.ratio p.value
 Left,Detail - Right,Detail    0.0382 0.0603   58.1  0.633  0.9209 
 Left,Detail - Left,Entire    -0.0782 0.0115 9217.0 -6.778  <.0001 
 Left,Detail - Right,Entire    0.1661 0.0603   58.1  2.753  0.0383 
 Right,Detail - Left,Entire   -0.1164 0.0603   58.1 -1.930  0.2268 
 Right,Detail - Right,Entire   0.1279 0.0115 9217.0 11.078  <.0001 
 Left,Entire - Right,Entire    0.2443 0.0603   58.1  4.050  0.0009 

Results are averaged over the levels of: Resolution 
P value adjustment: tukey method for comparing a family of 4 estimates 

2.1 Total Number of Fixations Hierarchical Steps —

Total Number of Fixations Nested Models

null <- nlme::gls(NFixAOIQRMNorm ~ 1, data = FixNumber)
AOI <- lm(NFixAOIQRMNorm ~ AOI, data = FixNumber)
AOIR <- lm(NFixAOIQRMNorm ~ AOI + Resolution, data = FixNumber)
AOIRM <- lm(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude, data = FixNumber)
AOIRM.intAOIM <- lm(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude, data = FixNumber)
AOIRM.intAOIM.rand.intAOIM <- lme4::lmer(NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixNumber)

Total Number of Fixations Nested Models AIC

MuMIn::AICc(null)

[1] 11932.01

MuMIn::AICc(AOI)

[1] 11703.74

MuMIn::AICc(AOIR)

[1] 11694.27

MuMIn::AICc(AOIRM)

[1] 11689.35

MuMIn::AICc(AOIRM.intAOIM)

[1] 11571.01

MuMIn::AICc(AOIRM.intAOIM.rand.intAOIM)

[1] 9274.262

Total Number of Fixations Nested Models BIC

BIC(null)

[1] 11946.28

BIC(AOI)

[1] 11725.14

BIC(AOIR)

[1] 11722.81

BIC(AOIRM)

[1] 11725.02

BIC(AOIRM.intAOIM)

[1] 11613.81

BIC(AOIRM.intAOIM.rand.intAOIM)

[1] 9324.198

Total Number of Fixations Nested Models R squared

MuMIn::r.squaredLR(null)

[1] -4.218847e-15
attr(,"adj.r.squared")
[1] -5.833118e-15

MuMIn::r.squaredGLMM(AOI)

            R2m        R2c
[1,] 0.02358046 0.02358046

MuMIn::r.squaredGLMM(AOIR)

            R2m        R2c
[1,] 0.02478398 0.02478398

MuMIn::r.squaredGLMM(AOIRM)

            R2m        R2c
[1,] 0.02550836 0.02550836

MuMIn::r.squaredGLMM(AOIRM.intAOIM)

            R2m        R2c
[1,] 0.03805875 0.03805875

MuMIn::r.squaredGLMM(AOIRM.intAOIM.rand.intAOIM)

            R2m       R2c
[1,] 0.03773117 0.2759189

Total Number of Fixations Nested Models Likelihood ratio test

lmtest::lrtest(null, AOI, AOIR, AOIRM, AOIRM.intAOIM, AOIRM.intAOIM.rand.intAOIM)

Likelihood ratio test

Model 1: NFixAOIQRMNorm ~ 1
Model 2: NFixAOIQRMNorm ~ AOI
Model 3: NFixAOIQRMNorm ~ AOI + Resolution
Model 4: NFixAOIQRMNorm ~ AOI + Resolution + Magnitude
Model 5: NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude
Model 6: NFixAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + 
    (1 | AOI:RecordingNameII)
  #Df  LogLik Df     Chisq Pr(>Chisq)    
1   2 -5964.0                            
2   3 -5848.9  1  230.2745  < 2.2e-16 ***
3   4 -5843.1  1   11.4690  0.0007077 ***
4   5 -5839.7  1    6.9211  0.0085185 ** 
5   6 -5779.5  1  120.3424  < 2.2e-16 ***
6   7 -4630.1  1 2298.7500  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

2.2 Total Number of Fixations Plots —————-

Interaction AOI * MAgnitude

Create Database with Total Number of Fixations fitted values (representing the effect on the scale of the linear predictor)

efNFix <- effects::effect("AOI:Magnitude", NFix)
summary(efNFix)

 AOI*Magnitude effect
       Magnitude
AOI        Detail    Entire
  Left  1.0332142 1.1114489
  Right 0.9950109 0.8671457

 Lower 95 Percent Confidence Limits
       Magnitude
AOI        Detail   Entire
  Left  0.9496052 1.027839
  Right 0.9114020 0.783536

 Upper 95 Percent Confidence Limits
       Magnitude
AOI       Detail    Entire
  Left  1.116823 1.1950585
  Right 1.078620 0.9507553

efNFix <- as.data.frame(efNFix)

Create Database with Total Number of Fixations means by Subject

require(dplyr)
NFixmeansSubj <- FixNumber %>% 
  group_by(AOI, Magnitude, RecordingNameII) %>% 
  summarise(RespM = mean(NFixAOIQRMNorm), Respse = sjstats::se(NFixAOIQRMNorm))
print(NFixmeansSubj)

2.2.0.1 Plot Total Number of Fixations fitted values

require(ggplot2)
NFixbar <- ggplot(efNFix, aes(AOI, fit, fill= Magnitude)) + 
  geom_point() + 
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4) + theme_bw(base_size=12)
NFixbar <- ggplot(efNFix, aes(AOI, fit, fill=Magnitude)) + geom_bar(stat="identity", position="dodge") +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4, position=position_dodge(width=0.9)) + 
  ylim(0,1.3) +
  labs(x= "AOI", y= "Normalized Total Fixations Number") +
  scale_fill_grey(name = "Magnitude") +
  ggtitle("AOI*Magnitude") +
  theme_classic() +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(NFixbar)

2.2.0.2 Plot Total Number of Fixations means by Subject

require(ggplot2)
NFixsplitSUBJ <- ggplot(NFixmeansSubj, aes(AOI, RespM, fill= Magnitude)) + 
    geom_bar(stat="identity", position="dodge") +
    facet_wrap(~ RecordingNameII) + 
    xlab("AOI") + 
    ylab("Normalized Total Fixations Number") +
    scale_fill_grey(name = "Magnitude", breaks=c("Full", "Detail"), labels=c("Entire", "Detail"))
print(NFixsplitSUBJ)

2.2.0.3 Boxplot Total Number of Fixations

NFixbox <- ggplot(data = FixNumber, aes(x = Magnitude, y = NFixAOIQRMNorm)) +
  facet_wrap(. ~ AOI) + 
  geom_jitter(alpha = 0.1, color = "black") + 
  geom_boxplot(alpha = 0) + 
  theme_light() + 
  labs(tag = "B", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized Total Fixations Number") +
  theme(plot.title = element_text(hjust = 0.5))
print(NFixbox)

2.2.0.4 Violinplot Total Number of Fixations

Function: plot split violins

GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
                             
                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })
geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}

Violinplot Total Number of Fixations

NFixsplitviolin <- ggplot(FixNumber, aes(x=Magnitude, y=NFixAOIQRMNorm, fill=AOI)) + 
  geom_split_violin() +
  scale_fill_grey() +
  theme_light() +
  labs(tag = "B", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized Total Fixations Number") +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(NFixsplitviolin)

3 Total Fixations Duration ———————–

Normalized Total Fixations Duration Database Import

FixDuration <- readr::read_delim("FixDuration.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)

Parsed with column specification:
cols(
  RecMedia = col_character(),
  RecordingNameII = col_character(),
  MediaName = col_character(),
  ImageStartTimestamp = col_double(),
  Cross = col_character(),
  Question = col_character(),
  Resolution = col_character(),
  Magnitude = col_character(),
  AOI = col_character(),
  GazeEventDuration = col_double(),
  DurAOIQRM = col_double(),
  DurAOIQRMNorm = col_double()
)

FixDuration <- subset(FixDuration, !(FixDuration$Question == "Baseline"))
FixDuration$AOI <- as.factor(FixDuration$AOI)
FixDuration$Magnitude <- as.factor(FixDuration$Magnitude)
levels(FixDuration$Magnitude) <- c("Detail", "Entire")
FixDuration$Resolution <- as.factor(FixDuration$Resolution)
print(FixDuration)

Dur Definitive Model

Dur <- lme4::lmer(DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixDuration)

Dur Overall Statistics

AIC(Dur)

[1] 11634.44

MuMIn::r.squaredGLMM(Dur)

            R2m       R2c
[1,] 0.03637073 0.2612974

Dur Summary & ANOVA

summary(Dur)

Linear mixed model fit by REML ['lmerMod']
Formula: DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude +      (1 | AOI:RecordingNameII)
   Data: FixDuration

REML criterion at convergence: 11620.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5739 -0.6694 -0.0575  0.6003  7.3576 

Random effects:
 Groups              Name        Variance Std.Dev.
 AOI:RecordingNameII (Intercept) 0.06071  0.2464  
 Residual                        0.19939  0.4465  
Number of obs: 9278, groups:  AOI:RecordingNameII, 58

Fixed effects:
                          Estimate Std. Error t value
(Intercept)               1.034596   0.046914  22.053
AOIRight                 -0.029767   0.066022  -0.451
ResolutionLow            -0.039760   0.009272  -4.288
MagnitudeEntire           0.143859   0.013112  10.972
AOIRight:MagnitudeEntire -0.241181   0.018543 -13.006

Correlation of Fixed Effects:
            (Intr) AOIRgh RsltnL MgntdE
AOIRight    -0.704                     
ResolutinLw -0.099  0.000              
MagnitdEntr -0.140  0.099  0.000       
AOIRght:MgE  0.099 -0.140  0.000 -0.707

aovDur <- aov(DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + Error(RecordingNameII/(Resolution + Magnitude)), data = FixDuration)
summary(aovDur)

Error: RecordingNameII
           Df Sum Sq Mean Sq F value Pr(>F)
Resolution  1  0.033  0.0328    0.04  0.843
Residuals  27 22.234  0.8235               

Error: RecordingNameII:Resolution
           Df Sum Sq Mean Sq F value   Pr(>F)    
Resolution  1  3.668   3.668  27.156 1.73e-05 ***
Magnitude   1  0.002   0.002   0.012    0.914    
Residuals  27  3.647   0.135                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: RecordingNameII:Magnitude
          Df Sum Sq Mean Sq F value  Pr(>F)   
Magnitude  1  1.256  1.2560   9.055 0.00549 **
Residuals 28  3.884  0.1387                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
                Df Sum Sq Mean Sq F value Pr(>F)    
AOI              1   52.4   52.45   203.9 <2e-16 ***
AOI:Magnitude    1   33.8   33.75   131.2 <2e-16 ***
Residuals     9189 2363.1    0.26                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Dur Robustlmm

require(lme4)
require(robustlmm)
robustDur <- robustlmm::rlmer(DurAOIQRMNorm ~  AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixDuration, rho.sigma.e = psi2propII(smoothPsi, k = 2.28), rho.sigma.b = psi2propII(smoothPsi, k = 2.28))

summary(robustDur)

Robust linear mixed model fit by DAStau 
Formula: DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude +      (1 | AOI:RecordingNameII) 
   Data: FixDuration 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.7653 -0.6558 -0.0338  0.6596  8.0822 

Random effects:
 Groups              Name        Variance Std.Dev.
 AOI:RecordingNameII (Intercept) 0.06492  0.2548  
 Residual                        0.18077  0.4252  
Number of obs: 9278, groups: AOI:RecordingNameII, 58

Fixed effects:
                          Estimate Std. Error t value
(Intercept)               1.012527   0.049572  20.426
AOIRight                 -0.020332   0.069812  -0.291
ResolutionLow            -0.042096   0.009054  -4.649
MagnitudeEntire           0.146255   0.012805  11.422
AOIRight:MagnitudeEntire -0.254199   0.018109 -14.037

Correlation of Fixed Effects:
            (Intr) AOIRgh RsltnL MgntdE
AOIRight    -0.704                     
ResolutinLw -0.091  0.000              
MagnitdEntr -0.129  0.092  0.000       
AOIRght:MgE  0.091 -0.130  0.000 -0.707

Robustness weights for the residuals: 
 7395 weights are ~= 1. The remaining 1883 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.166   0.670   0.822   0.786   0.939   0.999 

Robustness weights for the random effects: 
 47 weights are ~= 1. The remaining 11 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.516   0.695   0.894   0.816   0.946   0.980 

Rho functions used for fitting:
  Residuals:
    eff: smoothed Huber (k = 1.345, s = 10) 
    sig: smoothed Huber, Proposal II (k = 2.28, s = 10) 
  Random Effects, variance component 1 (AOI:RecordingNameII):
    eff: smoothed Huber (k = 1.345, s = 10) 
    vcp: smoothed Huber, Proposal II (k = 2.28, s = 10) 

Comparing lmer-rlmer models

campareDur <- robustlmm::compare(Dur, robustDur, show.rho.functions = FALSE)

Get the bootstrapped confidence intervals for the model parameters

ciDur <- lme4::confint.merMod(Dur, method = "boot", boot.type = "perc", nsim = 10000, oldNames = F)

ciDur

                                         2.5 %      97.5 %
sd_(Intercept)|AOI:RecordingNameII  0.26053283  0.38945797
sigma                               1.13069224  1.16372164
(Intercept)                         0.70949223  0.96761043
AOIRight                            0.03295627  0.39482166
AOILeft:MagnitudeEntire            -0.15037329 -0.01928073
AOIRight:MagnitudeEntire            0.09381820  0.22762033

Dur Post Hoc

emmeans::emm_options(pbkrtest.limit = 9278)
emmDur <- emmeans::emmeans(Dur, pairwise ~ AOI:Magnitude, adjust="tukey")

emmDur

$emmeans
 AOI   Magnitude emmean     SE   df lower.CL upper.CL
 Left  Detail     1.015 0.0467 58.3    0.921    1.108
 Right Detail     0.985 0.0467 58.3    0.892    1.078
 Left  Entire     1.159 0.0467 58.3    1.065    1.252
 Right Entire     0.888 0.0467 58.3    0.794    0.981

Results are averaged over the levels of: Resolution 
Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$contrasts
 contrast                    estimate     SE     df t.ratio p.value
 Left,Detail - Right,Detail    0.0298 0.0660   58.3   0.451 0.9692 
 Left,Detail - Left,Entire    -0.1439 0.0131 9217.0 -10.972 <.0001 
 Left,Detail - Right,Entire    0.1271 0.0660   58.3   1.925 0.2290 
 Right,Detail - Left,Entire   -0.1736 0.0660   58.3  -2.630 0.0518 
 Right,Detail - Right,Entire   0.0973 0.0131 9217.0   7.422 <.0001 
 Left,Entire - Right,Entire    0.2709 0.0660   58.3   4.104 0.0007 

Results are averaged over the levels of: Resolution 
P value adjustment: tukey method for comparing a family of 4 estimates 

3.1 Total Fixations Duration Hierarchical Steps —-

Total Fixations Duration Nested Models

null <- nlme::gls(DurAOIQRMNorm ~ 1, data = FixDuration)
AOI <- lm(DurAOIQRMNorm ~ AOI, data = FixDuration)
AOIR <- lm(DurAOIQRMNorm ~ AOI + Resolution, data = FixDuration)
AOIRM <- lm(DurAOIQRMNorm ~ AOI + Resolution + Magnitude, data = FixDuration)
AOIRM.intAOIM <- lm(DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude, data = FixDuration)
AOIRM.intAOIM.rand.intAOIM <- lme4::lmer(DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + (1|AOI:RecordingNameII), data = FixDuration)

Total Fixations Duration Nested Models AIC

MuMIn::AICc(null)

[1] 14116.18

MuMIn::AICc(AOI)

[1] 13911.57

MuMIn::AICc(AOIR)

[1] 13899.57

MuMIn::AICc(AOIRM)

[1] 13896.77

MuMIn::AICc(AOIRM.intAOIM)

[1] 13768.82

MuMIn::AICc(AOIRM.intAOIM.rand.intAOIM)

[1] 11634.46

Total Fixations Duration Nested Models BIC

BIC(null)

[1] 14130.45

BIC(AOI)

[1] 13932.98

BIC(AOIR)

[1] 13928.1

BIC(AOIRM)

[1] 13932.44

BIC(AOIRM.intAOIM)

[1] 13811.62

BIC(AOIRM.intAOIM.rand.intAOIM)

[1] 11684.39

Total Fixations Duration Nested Models R squared

MuMIn::r.squaredLR(null)

[1] 3.330669e-15
attr(,"adj.r.squared")
[1] 4.262919e-15

MuMIn::r.squaredGLMM(AOI)

            R2m        R2c
[1,] 0.02111166 0.02111166

MuMIn::r.squaredGLMM(AOIR)

          R2m      R2c
[1,] 0.022586 0.022586

MuMIn::r.squaredGLMM(AOIRM)

            R2m        R2c
[1,] 0.02308885 0.02308885

MuMIn::r.squaredGLMM(AOIRM.intAOIM)

            R2m        R2c
[1,] 0.03666878 0.03666878

MuMIn::r.squaredGLMM(AOIRM.intAOIM.rand.intAOIM)

            R2m       R2c
[1,] 0.03637073 0.2612974

Total Fixations Duration Nested Models Likelihood ratio test

lmtest::lrtest(null, AOI, AOIR, AOIRM, AOIRM.intAOIM, AOIRM.intAOIM.rand.intAOIM)

Likelihood ratio test

Model 1: DurAOIQRMNorm ~ 1
Model 2: DurAOIQRMNorm ~ AOI
Model 3: DurAOIQRMNorm ~ AOI + Resolution
Model 4: DurAOIQRMNorm ~ AOI + Resolution + Magnitude
Model 5: DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude
Model 6: DurAOIQRMNorm ~ AOI + Resolution + Magnitude + AOI:Magnitude + 
    (1 | AOI:RecordingNameII)
  #Df  LogLik Df     Chisq Pr(>Chisq)    
1   2 -7056.1                            
2   3 -6952.8  1  206.6075  < 2.2e-16 ***
3   4 -6945.8  1   14.0085   0.000182 ***
4   5 -6943.4  1    4.7986   0.028483 *  
5   6 -6878.4  1  129.9547  < 2.2e-16 ***
6   7 -5810.2  1 2136.3670  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

3.2 Total Fixations Duration Plots —————–

Create Database with Total Fixations Duration fitted values (representing the effect on the scale of the linear predictor) Interaction AOI * MAgnitude

efDur <- effects::effect("AOI:Magnitude", Dur)
summary(efDur)

 AOI*Magnitude effect
       Magnitude
AOI        Detail    Entire
  Left  1.0147200 1.1585791
  Right 0.9849531 0.8876315

 Lower 95 Percent Confidence Limits
       Magnitude
AOI        Detail    Entire
  Left  0.9232079 1.0670662
  Right 0.8934410 0.7961186

 Upper 95 Percent Confidence Limits
       Magnitude
AOI       Detail    Entire
  Left  1.106232 1.2500919
  Right 1.076465 0.9791443

efDur <- as.data.frame(efDur)

Plot Total Fixations Duration fitted values

require(ggplot2)
Durbar <- ggplot(efDur, aes(AOI, fit, fill= Magnitude)) + 
  geom_point() + 
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4) + theme_bw(base_size=12)
Durbar <- ggplot(efDur, aes(AOI, fit, fill=Magnitude)) + geom_bar(stat="identity", position="dodge") +
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4, position=position_dodge(width=0.9)) + 
  ylim(0,1.3) +
  labs(x= "AOI", y= "Normalized Total Fixations Duration") +
  scale_fill_grey(name = "Magnitude") +
  ggtitle("AOI*Magnitude") +
  theme_classic() +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(Durbar)

Create Database with Total Fixations Duration means by Subject

require(dplyr)
DurmeansSubj <- FixDuration %>% 
  group_by(AOI, Magnitude, RecordingNameII) %>% 
  summarise(RespM = mean(DurAOIQRMNorm), Respse = sjstats::se(DurAOIQRMNorm))
print(DurmeansSubj)

Plot Total Fixations Duration means by Subject

require(ggplot2)
DursplitSUBJ <- ggplot(DurmeansSubj, aes(AOI, RespM, fill= Magnitude)) + 
    geom_bar(stat="identity", position="dodge") +
    facet_wrap(~ RecordingNameII) + 
    xlab("AOI") + 
    ylab("Normalized Total Fixations Duration") +
    scale_fill_grey(name = "Magnitude", breaks=c("Full", "Detail"), labels=c("Entire", "Detail"))
print(DursplitSUBJ)

Boxplot Total Fixations Duration

Durbox <- ggplot(data = FixDuration, aes(x = Magnitude, y = DurAOIQRMNorm)) +
  facet_wrap(. ~ AOI) + 
  geom_jitter(alpha = 0.1, color = "black") + 
  geom_boxplot(alpha = 0) +
  theme_light() +
  labs(tag = "C", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized Total Fixations Duration") +
  theme(plot.title = element_text(hjust = 0.5))
print(Durbox)

Function: plot split violins

GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, 
                           draw_group = function(self, data, ..., draw_quantiles = NULL) {
                             data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth * (xmax - x))
                             grp <- data[1, "group"]
                             newdata <- plyr::arrange(transform(data, x = if (grp %% 2 == 1) xminv else xmaxv), if (grp %% 2 == 1) y else -y)
                             newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1, ])
                             newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
                             
                             if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
                               stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <=
                                                                         1))
                               quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
                               aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")), drop = FALSE]
                               aesthetics$alpha <- rep(1, nrow(quantiles))
                               both <- cbind(quantiles, aesthetics)
                               quantile_grob <- GeomPath$draw_panel(both, ...)
                               ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata, ...), quantile_grob))
                             }
                             else {
                               ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
                             }
                           })
geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity", ..., 
                              draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, 
                              show.legend = NA, inherit.aes = TRUE) {
  layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, 
        position = position, show.legend = show.legend, inherit.aes = inherit.aes, 
        params = list(trim = trim, scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}

Violinplot Total Fixations Duration

Dursplitviolin <- ggplot(FixDuration, aes(x=Magnitude, y=DurAOIQRMNorm, fill=AOI)) + 
  geom_split_violin() +
  scale_fill_grey() +
  theme_light() +
  labs(tag = "C", title = "AOI * Magnitude", x = "Magnitude", y = "Normalized Total Fixations Duration") +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
print(Dursplitviolin)

4 Behavioral ————————————-

Normalized Behavioral Scores Database Import

BehavioralD <- readr::read_delim("BehavioralD.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)

Parsed with column specification:
cols(
  RecMedia = col_character(),
  Question = col_character(),
  Resolution = col_character(),
  Magnitude = col_character(),
  RecordingNameII = col_character(),
  MediaName = col_character(),
  Resp = col_double()
)

BehavioralD <- subset(BehavioralD, !(BehavioralD$Question == "Baseline"))
BehavioralD <- subset(BehavioralD, !(is.na(BehavioralD$Resp)))
BehavioralD$Magnitude <- as.factor(BehavioralD$Magnitude)
levels(BehavioralD$Magnitude) <- c("Detail", "Entire")
print(BehavioralD)

Behavioral Cumulative Link Model for Ordinal Data (Likert)

Behavioral <- ordinal::clm(as.factor(Resp) ~ Question + Resolution + Magnitude + Question:Resolution + Question:Magnitude, data = BehavioralD, method = c("Newton"),  threshold="equidistant")

Behavioral Overall Statistics

AIC(Behavioral)

[1] 14657.83

Behavioral Summary & ANOVA

summary(Behavioral)

formula: as.factor(Resp) ~ Question + Resolution + Magnitude + Question:Resolution + Question:Magnitude
data:    BehavioralD

Coefficients:
                                   Estimate Std. Error z value Pr(>|z|)    
QuestionLike                       0.014546   0.142900   0.102 0.918924    
QuestionMovement                  -0.177560   0.145472  -1.221 0.222246    
QuestionProximity                 -0.013969   0.140602  -0.099 0.920860    
QuestionTouch                     -0.637380   0.144928  -4.398 1.09e-05 ***
ResolutionLow                     -0.944220   0.115561  -8.171 3.07e-16 ***
MagnitudeEntire                    0.055318   0.114824   0.482 0.629971    
QuestionLike:ResolutionLow         0.078474   0.162878   0.482 0.629950    
QuestionMovement:ResolutionLow     0.586653   0.165665   3.541 0.000398 ***
QuestionProximity:ResolutionLow    0.524013   0.160458   3.266 0.001092 ** 
QuestionTouch:ResolutionLow        0.335452   0.164202   2.043 0.041060 *  
QuestionLike:MagnitudeEntire      -0.111713   0.162785  -0.686 0.492550    
QuestionMovement:MagnitudeEntire  -0.003161   0.165319  -0.019 0.984743    
QuestionProximity:MagnitudeEntire -0.506254   0.160389  -3.156 0.001597 ** 
QuestionTouch:MagnitudeEntire     -0.307219   0.164063  -1.873 0.061128 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
            Estimate Std. Error z value
threshold.1 -1.76637    0.10486  -16.84
spacing      1.17395    0.01766   66.47

stats::anova(Behavioral, type="II")

Type II Analysis of Deviance Table with Wald chi-square tests

                    Df    Chisq Pr(>Chisq)    
Question             4 100.2918  < 2.2e-16 ***
Resolution           1 149.4088  < 2.2e-16 ***
Magnitude            1   6.2947  0.0121103 *  
Question:Resolution  4  20.3310  0.0004296 ***
Question:Magnitude   4  14.6347  0.0055221 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Behavioral Post Hoc

Interaction Question * Resolution

intQR <- emmeans::emmeans(Behavioral, pairwise ~  Question:Resolution, adjust="tukey")
# contrast                            estimate    SE  df z.ratio p.value
# Beautiful,Detail - Like,Detail      -0.05378 0.116 Inf -0.465  1.0000
# Beautiful,Detail - Movement,Detail  -0.11577 0.117 Inf -0.992  0.9928
# Beautiful,Detail - Proximity,Detail -0.24804 0.114 Inf -2.171  0.4770
# Beautiful,Detail - Touch,Detail      0.46965 0.117 Inf  4.009  0.0025 *
# Like,Detail - Movement,Detail       -0.06198 0.117 Inf -0.531  1.0000
# Like,Detail - Proximity,Detail      -0.19425 0.114 Inf -1.698  0.7971
# Like,Detail - Touch,Detail           0.52344 0.117 Inf  4.461  0.0003 *
# Movement,Detail - Proximity,Detail  -0.13227 0.115 Inf -1.146  0.9799
# Movement,Detail - Touch,Detail       0.58542 0.118 Inf  4.949  <.0001 *
# Proximity,Detail - Touch,Detail      0.71769 0.116 Inf  6.185  <.0001 *
# 
# Beautiful,Detail - Beautiful,Full   -0.05532 0.115 Inf -0.482  1.0000
# Like,Detail - Like,Full              0.05639 0.115 Inf  0.489  1.0000
# Movement,Detail - Movement,Full     -0.05216 0.119 Inf -0.439  1.0000
# Proximity,Detail - Proximity,Full    0.45094 0.112 Inf  4.027  0.0023 *
# Touch,Detail - Touch,Full            0.25190 0.117 Inf  2.150  0.4920
# 
# Beautiful,Full - Like,Full           0.05793 0.115 Inf  0.506  1.0000
# Beautiful,Full - Movement,Full      -0.11261 0.117 Inf -0.960  0.9943
# Beautiful,Full - Proximity,Full      0.25822 0.113 Inf  2.295  0.3923
# Beautiful,Full - Touch,Full          0.77687 0.115 Inf  6.730  <.0001 *
# Like,Full - Movement,Full           -0.17054 0.118 Inf -1.450  0.9113
# Like,Full - Proximity,Full           0.20029 0.113 Inf  1.775  0.7515
# Like,Full - Touch,Full               0.71894 0.116 Inf  6.212  <.0001 *
# Movement,Full - Proximity,Full       0.37082 0.116 Inf  3.207  0.0438
# Movement,Full - Touch,Full           0.88948 0.118 Inf  7.508  <.0001 *
# Proximity,Full - Touch,Full          0.51866 0.113 Inf  4.573  0.0002 *

Interaction Question * Magnitude

intQM <- emmeans::emmeans(Behavioral, pairwise ~  Question:Magnitude, adjust="tukey")
# contrast                        estimate    SE  df z.ratio p.value
# Beautiful,High - Like,High        0.0413 0.117 Inf  0.353  1.0000
# Beautiful,High - Movement,High    0.1791 0.120 Inf  1.495  0.8946
# Beautiful,High - Proximity,High   0.2671 0.115 Inf  2.318  0.3770
# Beautiful,High - Touch,High       0.7910 0.118 Inf  6.682  <.0001 *
# Like,High - Movement,High         0.1378 0.120 Inf  1.152  0.9791
# Like,High - Proximity,High        0.2258 0.115 Inf  1.963  0.6252
# Like,High - Touch,High            0.7497 0.118 Inf  6.347  <.0001 *
# Movement,High - Proximity,High    0.0880 0.118 Inf  0.747  0.9992
# Movement,High - Touch,High        0.6118 0.121 Inf  5.073  <.0001 *
# Proximity,High - Touch,High       0.5239 0.116 Inf  4.517  0.0003 *
# 
# Beautiful,High - Beautiful,Low    0.9442 0.116 Inf  8.171  <.0001 *
# Like,High - Like,Low              0.8657 0.116 Inf  7.470  <.0001 *
# Movement,High - Movement,Low      0.3576 0.119 Inf  3.001  0.0802
# Proximity,High - Proximity,Low    0.4202 0.112 Inf  3.754  0.0067 *
# Touch,High - Touch,Low            0.6088 0.117 Inf  5.188  <.0001 *
# 
# Beautiful,Low - Like,Low         -0.0372 0.113 Inf -0.328  1.0000
# Beautiful,Low - Movement,Low     -0.4075 0.114 Inf -3.566  0.0134 *
# Beautiful,Low - Proximity,Low    -0.2569 0.112 Inf -2.302  0.3872
# Beautiful,Low - Touch,Low         0.4555 0.114 Inf  3.986  0.0027 *
# Like,Low - Movement,Low          -0.3703 0.115 Inf -3.222  0.0419
# Like,Low - Proximity,Low         -0.2198 0.112 Inf -1.957  0.6292
# Like,Low - Touch,Low              0.4927 0.115 Inf  4.285  0.0008 *
# Movement,Low - Proximity,Low      0.1506 0.113 Inf  1.330  0.9470
# Movement,Low - Touch,Low          0.8631 0.116 Inf  7.429  <.0001 *
# Proximity,Low - Touch,Low         0.7125 0.113 Inf  6.280  <.0001 *

4.1 Behavioral Hierarchical Steps ——————

BehavioralNested Models

null <- ordinal::clm(as.factor(Resp) ~ 1, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
Q <- ordinal::clm(as.factor(Resp)  ~ Question, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
QR <- ordinal::clm(as.factor(Resp)  ~ Question + Resolution, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
QRM <- ordinal::clm(as.factor(Resp)  ~ Question + Resolution + Magnitude, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
QRM.intQR <- ordinal::clm(as.factor(Resp)  ~ Question + Resolution + Magnitude + Question:Resolution, data = BehavioralD, method = c("Newton"),  threshold="equidistant")
QRM.intQR.intQM <- ordinal::clm(as.factor(Resp)  ~ Question + Resolution + Magnitude + Question:Resolution +  Question:Magnitude, data = BehavioralD, method = c("Newton"),  threshold="equidistant")

Behavioral Nested Models AIC

MuMIn::AICc(null)

[1] 14920.79

MuMIn::AICc(Q)

[1] 14830.54

MuMIn::AICc(QR)

[1] 14681.33

MuMIn::AICc(QRM) 

[1] 14676.55

MuMIn::AICc(QRM.intQR)

[1] 14664.54

MuMIn::AICc(QRM.intQR.intQM)

[1] 14657.94

Behavioral Nested Models Likelihood ratio test

lmtest::lrtest(null, Q, QR, QRM, QRM.intQR, QRM.intQR.intQM)

Likelihood ratio test

Model 1: as.factor(Resp) ~ 1
Model 2: as.factor(Resp) ~ Question
Model 3: as.factor(Resp) ~ Question + Resolution
Model 4: as.factor(Resp) ~ Question + Resolution + Magnitude
Model 5: as.factor(Resp) ~ Question + Resolution + Magnitude + Question:Resolution
Model 6: as.factor(Resp) ~ Question + Resolution + Magnitude + Question:Resolution + 
    Question:Magnitude
  #Df  LogLik Df    Chisq Pr(>Chisq)    
1   2 -7458.4                           
2   6 -7409.3  4  98.2687  < 2.2e-16 ***
3   7 -7333.7  1 151.2154  < 2.2e-16 ***
4   8 -7330.3  1   6.7828   0.009204 ** 
5  12 -7320.2  4  20.0462   0.000489 ***
6  16 -7312.9  4  14.6474   0.005491 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4.2 Behavioral Plots ——————————-

4.2.0.1 Diverging stacked bar chart Question * Resolution

QRLikert <- readr::read_delim("QRLikert.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)

Parsed with column specification:
cols(
  Beautiful.High = col_double(),
  Like.High = col_double(),
  Movement.High = col_double(),
  Proximity.High = col_double(),
  Touch.High = col_double(),
  Beautiful.Low = col_double(),
  Like.Low = col_double(),
  Movement.Low = col_double(),
  Proximity.Low = col_double(),
  Touch.Low = col_double()
)

require(sjPlot)

Carico il pacchetto richiesto: sjPlot
Package version inconsistency detected.
TMB was built with Matrix version 1.2.17
Current Matrix version is 1.2.15
Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package

require(sjmisc)

Carico il pacchetto richiesto: sjmisc

levels_5 <- c("1", "2", "3", "4", "5")
itemsQR <- c("Beautiful-High","Like-High","Movement-High","Proximity-High","Touch-High",
              "Beautiful-Low","Like-Low","Movement-Low","Proximity-Low","Touch-Low")
plot_likert(QRLikert,
            expand.grid = FALSE,
            axis.labels = itemsQR,
            legend.labels = levels_5,
            cat.neutral = 3,
            title = "Question * Resolution",
            legend.title = "Likert Scores",
            show.n = FALSE,
            coord.flip = T) + 
  font_size(labels.y = 13, title = 15, base.theme = legend_style(pos = "bottom", justify = "center", base.theme = theme_sjplot(base_size = 10)))

4.2.0.2 Diverging stacked bar chart Question * Magnitude

QMLikert <- readr::read_delim("QMLikert.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)

Parsed with column specification:
cols(
  Beautiful.Detail = col_double(),
  Like.Detail = col_double(),
  Movement.Detail = col_double(),
  Proximity.Detail = col_double(),
  Touch.Detail = col_double(),
  Beautiful.Entire = col_double(),
  Like.Entire = col_double(),
  Movement.Entire = col_double(),
  Proximity.Entire = col_double(),
  Touch.Entire = col_double()
)

require(sjPlot)
require(sjmisc)
levels_5 <- c("1", "2", "3", "4", "5")
itemsQM <- c("Beautiful-Detail","Like-Detail","Movement-Detail","Proximity-Detail","Touch-Detail",
             "Beautiful-Entire","Like-Entire","Movement-Entire","Proximity-Entire","Touch-Entire")
plot_likert(QMLikert,
            expand.grid = FALSE,
            axis.labels = itemsQM,
            legend.labels = levels_5,
            cat.neutral = 3,
            title = "Question * Magnitude",
            legend.title = "Likert Scores",
            show.n = FALSE,
            coord.flip = T) + 
  font_size(labels.y = 13, title = 15, base.theme = legend_style(pos = "bottom", justify = "center", base.theme = theme_sjplot(base_size = 10)))

4.3 Behavioral Corr ——————————–

Normalized Behavioral Scores means by Subject Dataset Import

Correlations <- readr::read_delim("Correlations.txt", "\t", escape_double = FALSE,  trim_ws = TRUE)

Parsed with column specification:
cols(
  Rec = col_double(),
  IRI = col_double(),
  IRIF.z = col_double(),
  BeautifulMeansNorm = col_double(),
  LikeMeansNorm = col_double(),
  MovementMeansNorm = col_double(),
  ProximityMeansNorm = col_double(),
  TouchMeansNorm = col_double()
)

print(Correlations)

Correlation Among Questions (Bonferroni’s Correction -> 0.05/10 = 0.005)

BLc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$LikeMeansNorm, method = c("spearman"))

Impossibile calcolare p-value esatti in presenza di ties

ifelse(BLc$p.value < 0.005, print(BLc), NA)

    Spearman's rank correlation rho

data:  Correlations$BeautifulMeansNorm and Correlations$LikeMeansNorm
S = 262.53, p-value = 9.189e-15
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.9415953 

[[1]]
       S 
262.5292 

BMc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$MovementMeansNorm, method = c("spearman"))
ifelse(BMc$p.value < 0.005, print(BMc), NA)

[1] NA

BPc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$ProximityMeansNorm, method = c("spearman"))
ifelse(BPc$p.value < 0.005, print(BPc), NA)

[1] NA

BTc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$TouchMeansNorm, method = c("spearman"))
ifelse(BTc$p.value < 0.005, print(BTc), NA)

    Spearman's rank correlation rho

data:  Correlations$BeautifulMeansNorm and Correlations$TouchMeansNorm
S = 1528, p-value = 0.0001062
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.6600667 

[[1]]
   S 
1528 

LMc <- cor.test(Correlations$LikeMeansNorm, Correlations$MovementMeansNorm, method = c("spearman"))

Impossibile calcolare p-value esatti in presenza di ties

ifelse(LMc$p.value < 0.005, print(LMc), NA)

[1] NA

LPc <- cor.test(Correlations$LikeMeansNorm, Correlations$ProximityMeansNorm, method = c("spearman"))

Impossibile calcolare p-value esatti in presenza di ties

ifelse(LPc$p.value < 0.005, print(LPc), NA)

[1] NA

LTc <- cor.test(Correlations$LikeMeansNorm, Correlations$TouchMeansNorm, method = c("spearman"))

Impossibile calcolare p-value esatti in presenza di ties

ifelse(LTc$p.value < 0.005, print(LTc), NA)

    Spearman's rank correlation rho

data:  Correlations$LikeMeansNorm and Correlations$TouchMeansNorm
S = 1285.6, p-value = 9.405e-06
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.7139838 

[[1]]
       S 
1285.643 

MPc <- cor.test(Correlations$MovementMeansNorm, Correlations$ProximityMeansNorm, method = c("spearman"))
ifelse(MPc$p.value < 0.005, print(MPc), NA)

[1] NA

MTc <- cor.test(Correlations$MovementMeansNorm, Correlations$TouchMeansNorm, method = c("spearman"))
ifelse(MTc$p.value < 0.005, print(MTc), NA)

[1] NA

PTc <- cor.test(Correlations$ProximityMeansNorm, Correlations$TouchMeansNorm, method = c("spearman"))
ifelse(PTc$p.value < 0.005, print(PTc), NA)

[1] NA

Correlationtion Questions * IRI Fantasy Scale (Bonferroni’s Correction -> 0.05/5 = 0.01)

BIc <- cor.test(Correlations$BeautifulMeansNorm, Correlations$IRIF.z, method = c("spearman"))

Impossibile calcolare p-value esatti in presenza di ties

ifelse(BIc$p.value < 0.01, print(BIc), NA)

[1] NA

MIc <- cor.test(Correlations$MovementMeansNorm, Correlations$IRIF.z, method = c("spearman"))

Impossibile calcolare p-value esatti in presenza di ties

ifelse(MIc$p.value < 0.01, print(MIc), NA)

[1] NA

PIc <- cor.test(Correlations$ProximityMeansNorm, Correlations$IRIF.z, method = c("spearman"))

Impossibile calcolare p-value esatti in presenza di ties

ifelse(PIc$p.value < 0.01, print(PIc), NA)

[1] NA

TIc <- cor.test(Correlations$TouchMeansNorm, Correlations$IRIF.z, method = c("spearman"))

Impossibile calcolare p-value esatti in presenza di ties

ifelse(TIc$p.value < 0.01, print(TIc), NA)

[1] NA

LIc <- cor.test(Correlations$LikeMeansNorm, Correlations$IRIF.z, method = c("spearman"))

Impossibile calcolare p-value esatti in presenza di ties

ifelse(LIc$p.value < 0.01, print(LIc), NA)

[1] NA

4.3.1 Plot Stat Significative Correlations

4.3.1.1 Like ~ Beautiful

require(ggplot2)
LB <- ggplot(Correlations, aes(x=LikeMeansNorm, y=BeautifulMeansNorm)) +
  geom_point(size=2) +
  geom_smooth(method=lm, se=T, color="black") +
  ylim(0.3, 1.5) +
  xlim(0.5, 1.4) +
  labs(x= "Like (mean score)", y= "Beautiful (mean score)") +
  ggtitle("") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))
print(LB)

4.3.1.2 Like ~ Touch

require(ggplot2)
LT <- ggplot(Correlations, aes(x=LikeMeansNorm, y=TouchMeansNorm)) +
  geom_point(size=2) +
  geom_smooth(method=lm, se=T, color="black") +
  ylim(0.3, 1.5) +
  xlim(0.5, 1.4) +
  labs(x= "Like (mean score)", y= "Touch (mean score)") +
  ggtitle("") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))
print(LT)

4.3.1.3 Beautiful ~ Touch

require(ggplot2)
BT <- ggplot(Correlations, aes(x=BeautifulMeansNorm, y=TouchMeansNorm)) +
  geom_point(size=2) +
  geom_smooth(method=lm, se=T, color="black") +
  ylim(0.3, 1.5) +
  xlim(0.5, 1.4) +
  labs(x= "Beautiful (mean score)", y= "Touch (mean score)") +
  ggtitle("") +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))
print(BT)

5 Citations ————————————–

Print References

General Packages

citation(package ='readr')

To cite package ‘readr’ in publications use:

  Hadley Wickham, Jim Hester and Romain Francois (2018). readr: Read Rectangular Text Data. R package version 1.3.1.
  https://CRAN.R-project.org/package=readr

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {readr: Read Rectangular Text Data},
    author = {Hadley Wickham and Jim Hester and Romain Francois},
    year = {2018},
    note = {R package version 1.3.1},
    url = {https://CRAN.R-project.org/package=readr},
  }

citation(package ='sjstats')

Lüdecke D (2019). _sjstats: Statistical Functions for Regression Models (Version 0.17.6)_. doi: 10.5281/zenodo.1284472 (URL:
http://doi.org/10.5281/zenodo.1284472), <URL: https://CRAN.R-project.org/package=sjstats>.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {sjstats: Statistical Functions for Regression Models (Version 0.17.6)},
    author = {Daniel Lüdecke},
    year = {2019},
    url = {https://CRAN.R-project.org/package=sjstats},
    doi = {10.5281/zenodo.1284472},
  }

citation(package ='MuMIn')

To cite package ‘MuMIn’ in publications use:

  Kamil Barton (2018). MuMIn: Multi-Model Inference. R package version 1.42.1. https://CRAN.R-project.org/package=MuMIn

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {MuMIn: Multi-Model Inference},
    author = {Kamil Barton},
    year = {2018},
    note = {R package version 1.42.1},
    url = {https://CRAN.R-project.org/package=MuMIn},
  }

ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may need manual editing, see
‘help("citation")’.

citation(package ='apaTables')

To cite package ‘apaTables’ in publications use:

  David Stanley (2018). apaTables: Create American Psychological Association (APA) Style Tables. R package version 2.0.5.
  https://CRAN.R-project.org/package=apaTables

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {apaTables: Create American Psychological Association (APA) Style Tables},
    author = {David Stanley},
    year = {2018},
    note = {R package version 2.0.5},
    url = {https://CRAN.R-project.org/package=apaTables},
  }

citation(package ='dplyr')

To cite package ‘dplyr’ in publications use:

  Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2019). dplyr: A Grammar of Data Manipulation. R package version
  0.8.3. https://CRAN.R-project.org/package=dplyr

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {dplyr: A Grammar of Data Manipulation},
    author = {Hadley Wickham and Romain François and Lionel Henry and Kirill Müller},
    year = {2019},
    note = {R package version 0.8.3},
    url = {https://CRAN.R-project.org/package=dplyr},
  }

Packages to cite for Statistics

citation(package ='lme4')

To cite lme4 in publications use:

  Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of
  Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.

A BibTeX entry for LaTeX users is

  @Article{,
    title = {Fitting Linear Mixed-Effects Models Using {lme4}},
    author = {Douglas Bates and Martin M{\"a}chler and Ben Bolker and Steve Walker},
    journal = {Journal of Statistical Software},
    year = {2015},
    volume = {67},
    number = {1},
    pages = {1--48},
    doi = {10.18637/jss.v067.i01},
  }

citation(package ='ez')

To cite package ‘ez’ in publications use:

  Michael A. Lawrence (2016). ez: Easy Analysis and Visualization of Factorial Experiments. R package version 4.4-0.
  https://CRAN.R-project.org/package=ez

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {ez: Easy Analysis and Visualization of Factorial Experiments},
    author = {Michael A. Lawrence},
    year = {2016},
    note = {R package version 4.4-0},
    url = {https://CRAN.R-project.org/package=ez},
  }

ATTENTION: This citation information has been auto-generated from the package DESCRIPTION file and may need manual editing, see
‘help("citation")’.

citation(package ='emmeans')

To cite package ‘emmeans’ in publications use:

  Russell Lenth (2019). emmeans: Estimated Marginal Means, aka Least-Squares Means. R package version 1.3.3.
  https://CRAN.R-project.org/package=emmeans

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {emmeans: Estimated Marginal Means, aka Least-Squares Means},
    author = {Russell Lenth},
    year = {2019},
    note = {R package version 1.3.3},
    url = {https://CRAN.R-project.org/package=emmeans},
  }

citation(package = 'robustlmm')

To cite robustlmm in publications use:

  Manuel Koller (2016). robustlmm: An R Package for Robust Estimation of Linear Mixed-Effects Models. Journal of Statistical Software,
  75(6), 1-24. doi:10.18637/jss.v075.i06

A BibTeX entry for LaTeX users is

  @Article{,
    title = {{robustlmm}: An {R} Package for Robust Estimation of Linear Mixed-Effects Models},
    author = {Manuel Koller},
    journal = {Journal of Statistical Software},
    year = {2016},
    volume = {75},
    number = {6},
    pages = {1--24},
    doi = {10.18637/jss.v075.i06},
  }

citation(package = 'ordinal')

To cite the ordinal-package in publications use:

  Christensen, R. H. B. (2019). ordinal - Regression Models for Ordinal Data. R package version 2019.4-25.
  http://www.cran.r-project.org/package=ordinal/.

A BibTeX entry for LaTeX users is

  @Misc{,
    title = {ordinal---Regression Models for Ordinal Data },
    author = {R. H. B. Christensen},
    year = {2019},
    note = {R package version 2019.4-25. http://www.cran.r-project.org/package=ordinal/},
  }

Packages to cite for Plotting

citation(package = 'effects')

To cite effects in publications use:

  John Fox and Sanford Weisberg (2019). An R Companion to Applied Regression, 3rd Edition. Thousand Oaks, CA
  <http://tinyurl.com/carbook>

To cite effects in publications use:

Fox J, Weisberg S (2018). “Visualizing Fit and Lack of Fit in Complex Regression Models with Predictor Effect Plots and Partial
Residuals.” _Journal of Statistical Software_, *87*(9), 1-27. doi: 10.18637/jss.v087.i09 (URL: http://doi.org/10.18637/jss.v087.i09).

For generalized linear models also cite:

  John Fox (2003). Effect Displays in R for Generalised Linear Models. Journal of Statistical Software, 8(15), 1-27. URL
  http://www.jstatsoft.org/v08/i15/.

For usage in multinomial and proportional-odds logit models also cite:

  John Fox, Jangman Hong (2009). Effect Displays in R for Multinomial and Proportional-Odds Logit Models: Extensions to the effects
  Package. Journal of Statistical Software, 32(1), 1-24. URL http://www.jstatsoft.org/v32/i01/.

To see these entries in BibTeX format, use 'print(<citation>, bibtex=TRUE)', 'toBibtex(.)', or set 'options(citation.bibtex.max=999)'.

citation(package ='ggplot2')

To cite ggplot2 in publications, please use:

  H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.

A BibTeX entry for LaTeX users is

  @Book{,
    author = {Hadley Wickham},
    title = {ggplot2: Elegant Graphics for Data Analysis},
    publisher = {Springer-Verlag New York},
    year = {2016},
    isbn = {978-3-319-24277-4},
    url = {https://ggplot2.tidyverse.org},
  }